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ABSTRACT 

. We develop a new method to determine the linear mass power spectrum using the 

£vq ■ mass function of galaxy clusters. We obtain the rms mass fluctuation a{M) using the 

expression for the mass function in the Press & Schechter (1974), Sheth, Mo & Tormen 
(2001) & Jenkins et al. (2001) formalisms. We apply different techniques to recover the 
adimensional power spectrum A 2 (A:) from a(M) namely the k e ff approximation, the 
singular value decomposition and the linear regularization method. The application 
of these techniques to the rCDM and ACDM GIF simulations shows a high efficiency 
in recovering the theoretical power spectrum over a wide range of scales. We compare 
our results with those derived from the power spectrum of the spatial distribution of 
the same sample of clusters in the simulations obtained by application of the classical 
| Feldman, Kaiser & Peacock (1994), FKP, method. We find that the mass function 

i based method presented here can provide a very accurate estimate of the linear power 

' spectrum, particularly for low values of k. This estimate is comparable, or even better 

behaved, than the FKP solution. 
Q ■ The principal advantage of our method is that it allows the determination of the 

S-^ | linear mass power spectrum using the joint information of objects of a wide range of 

c/3 , masses without dealing with specific assumptions on the bias relative to the underlying 

mass distribution. 
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1 INTRODUCTION the results obtained are affected by the incompleteness of 

the survey, its geometry, the redshift space distortions and 

Cosmological models predict the large scale matter distribu- the biag factor (which can be scale dependent) Moreover the 

tion evolved from primordial fluctuations in the density of distribution of galaxies or clusters of galaxies is affected by 

the universe. Unfortunately this distribution is not directly nondinear evolution effects, which makes it even more diffi- 

accessible through observations because we can only observe cmt tQ carry Qut a direct comparison between obse rvational 

the distributions of objects like galaxies or clusters of galax- , ,. , , , ,. .. 

J to b results and theoretical model predictions, 
ies that do not, necessarily, trace the mass distribution in a 

simple way. An alternative technique was presented by Gaztanaga & 

The power spectrum of the galaxy distribution has Baugh (1998), and has been applied in different ways by sev- 

proven to be the most popular statistic used to characterize eral other authors ( Dodelson & Gastanaga, 2001, Eisenstein 

the matter distribution in the universe. This is due to its & Zaldarriaga, 2001 and references therein). In it the power 

tight relation to physical properties of the growth of struc- spectrum of the galaxy distribution is derived from the an- 

tures throughout the history of the Universe. Also, the many gular correlation function inverting Limber's equation. As 

ways and different kind of data that can be used to measure this procedure does not involve any Fourier transforms nor 

the power spectrum has allowed independent results from the use of redshift data, the results obtained in this way 

redshift and angular galaxy surveys, cluster surveys, and are not affected by a convolution with the survey window 

CMB fluctuations studies. Each of these approaches have function, nor by redshift space distortions. The principal 

different drawbacks, though. For instance, the results from uncertainty in this approach comes from the fact that one 

large redshift surveys, are subject to several problems since of the ingredients of the technique is the redshift distribu- 



© 0000 RAS 



2 A. G. Sanchez, N. D. Padilla & D. G. Lambas 



tion of galaxies, which may not be accurately known, and 
moreover the results are still affected by the bias factor and 
non- linear evolution. 

The aim of this work is to develop a new method for 
obtaining the power spectrum, free of the problems of the 
methods previously described. In order to do this, we use the 
formalisms of Press & Schechter (1974, hereafter PS) and 
Sheth, Mo & Tormen (2001, hereafter SMT) where halos 
form in the higher peaks of the evolving primordial density 
fluctuations which are related to the linear power spectrum. 
This indicates that studying the distribution of halo masses, 
one can obtain information about the mass distribution, and 
therefore, recover the mass power spectrum. More precisely, 
PS's like recipes are sets of equations that allow the deter- 
mination of the mass function of dark matter halos from 
the mass power spectrum. If the mass function can be de- 
termined observationally, we can use these equations in the 
inverse way to derive the linear power spectrum. 

In this paper we present a new method for determining 
the linear power spectrum from the mass function of clusters 
of galaxies based on the PS, SMT and Jenkins et al. (2001, 
hereafter J01) prescriptions for the mass function of systems 
expected in hierarchical clustering scenarios. 

The outline of the paper is as follows, in section §2 we 
briefly review the different formalisms for obtaining a mass 
function used in this work. In §3 we give the general outline 
of our method and a simple test in an idealized case. In §4 
we describe a test of our method by its application to the 
GIF simulations ( Kauffman et al., 1999). In §4 we show 
a comparison with the results obtained using the standard 
technique of Feldman, Kaiser & Peacock (1994, hereafter 
FKP) for measuring the power spectrum, and finally, in §5 
we present a short discussion and the main conclusions. 



2 THE MASS FUNCTION 

An important tool with which to characterize the density 
field is the rms density fluctuation, smoothed on some co- 
moving scale R. For the case of a top-hat window, this scale 
is related to a mass M = ^-pR 3 , where p is the mean den- 
sity of the universe (p — p c O,). In turn, this quantity and 
the power spectrum P(k) = \Sk\ 2 are related by 



cr(R) 



A 2 (k)W 2 (kR) dlnfc 



(1) 



where A 2 (fc) = "^^i^) * s the dimensionless power spec- 
trum and W(kR) is the Fourier transform of the top-hat 
window function 



W(kR) = —MkR) 



,(sin(fei?) - kRcos(kR)) 

' (kkj* ' 



where ji(x) is the spherical Bessel function of order 1. 

PS developed a simple formalism to obtain an analytical 
expression for the differential mass function of dark matter 
halos based on the spherical collapse model. In this model, 
the abundance of haloes with masses M — > M + dM is given 
by: 



n{M)AM 



M 2 a(M) 



dln(a(M)) 



dln(M) 



(2) 
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Figure 1. The mass functions predicted by different formalisms 
for a -rCDM model. The curve obtained from the SMT formal- 
ism gives a much better description of the results from numerical 
simulations (shown by the J01 general fit) than the PS recipe. 



X exp 



2<t(M) 2 



dM, 



where 5 C is the density threshold corresponding to collapse 
according to the spherical collapse model. For an Q = 1 
universe 5 C = 1.686 and its value has little dependence on 
cosmology (Eke et al., 1996). This function is plotted in 
Figure 1 for a rCDM model. 

The vague analytical justification in the derivation of 
expression (El) has lead to different criticisms of the PS for- 
malism ( e.g. Peacock & Heavens, 1989), but a later and 
more formal derivation due to Bond et al. ( 1991) based on 
a particular type of filtering, and the set excursion theory, 
yields the same expression for the mass function according 
to the spherical collapse model, placing the PS formalism 
on a more solid theoretical basis. In addition to the original 
tests presented in PS, equation (Q) has been checked against 
results from a variety of numerical simulations (Efstathiou 
et al., 1988, Colberg & Couchman, 1989) showing a good 
general agreement. 

Nevertheless, recent works (Tormen, 1998; Sheth & Tor- 
men, 1999, Governato et al. 1999, J01) detected discrepan- 
cies between the PS predictions and the halo abundances 
obtained in much better numerical simulations. The PS for- 
mula underestimates the abundances of massive halos and 
overestimates the low mass end of the mass function (al- 
ready seen in the earlier numerical simulations results). A 
better description of the numerical results is given by a cor- 
rection to the PS mass function proposed originally by Sheth 
& Tormen (1999) 



n{M)AM = 



2aA 2 p 



7T M 2 <j{M) 



1 + 



\fa8 c J 



2P 



(3) 
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din (a) 



dln(M) 



( aS " \ 



AM 



(4) 



where A = 0.322, a = 0.707 and p = 0.3. This function is 
also plotted in Figure 1. SMT showed that this correction 
can be understood as a result of the incorporation of the 
ellipsoidal rather than spherical collapse model dynamics in 
a formalism similar to the one used by Bond et al. (1991). 
We can define a function 



,, , M dn cum {M,z) 



ciln<7~ 



(•») 



where n cum (M, z) is the cumulative mass function. From 
and (^) we can calculate f(a, z) for PS and SMT formalisms 
where 



fa; PS) = 

f(a;SMT) = 



2 8 C 



■ exp 



1 + 



A. 

2a 2 



a 

aS 2 



■ exp 



2a 2 



(6) 



(7) 



That is, the PS and SMT formalisms predict that the func- 
tion f(a, z) is independent of redshift and the cosmological 
model when it is expressed in terms of cr. 

From the analysis of a set of simulations carried out by 
the Virgo Consortium, J01 found that the function /(cr) can 
be described with an accuracy better than 20% for different 
redshifts and initial conditions by the general fit 



f(a) = A exp (- In cr -1 + B\ 



(8) 



where A = 0.315, B = 0.61 and e = 3.8, for halos identified 
with the FOF algorithm using the same linking length b = 
0.2 for every value of Q. 

The expression for n(M) inferred from (B) and (^) 



n(M,z)dM 



Ap(z) ri Infer' 1 ) 
M 2 dln(M) 



-(|m (CT -i+sn 



AM 



(9) 



This differential mass function is also shown in Figure 1, 
where it can be seen that SMT formalism is in a much better 
agreement with the general fit of J01 (that is, the results of 
numerical simulations), than PS formalism. 



3 THE METHOD 

In this section we analyse the problem of recovering the 
power spectrum A 2 (fc) using equation (|l|) and the expres- 
sions for the mass function in the formalisms just described. 
The steps to be followed in obtaining the linear mass power 
spectrum are: 

(i) Obtain the variance a(M) from the mass function. In 
order to do this we need to rewrite the expresions relating 
these statistics. 

(ii) Invert the power spectrum from the variance. This 
involves the application and test of different inversion algo- 
rithms. 



3.1 Obtaining a(M) from n(M) 

Assuming the mass function of dark matter haloes obeys one 
of the three formalisms described in the last section, we can 



re-interpret equations (^|), (||) and Q as differential equa- 
tions and use them in the determination of er(M) provided 
n(M) is known. In the Appendix we describe the solution 
of the equations obtained for each formalism. 

In the case of PS formalism the final solution is 

<r(M) = — T , (10) 



v^erT 



erf( 



where erf(:r) is the error function and G(M) is given by 



G(M) = 



A/8 



Mn(M)AM. 



a(M) 



For the SMT formalism we get 

Sc 



<5 C 



G(M) 
pA 



where we keep the definition of G(M) and 
$Or)=erf(x) + -=r(--pjp[- 



p, x 



(11) 



(12) 



where p is a parameter in the SMT formalism, and F (x) and 
P (a, x) are the gamma and incomplete gamma functions. 
Finally, for the J01 prescription, 



a(M) = exp 



*(ln (jg 



G(M) 



Ap 



(13) 



where the function *B(S) is given by 
i-s 

*(S) = / e- ,x ^ r Ax. 



In every case the values of Q and the amplitude of density 
fluctuations measured in a sphere of 8/i." 1 Mpc, erg, are as- 
sumed to be known. 

If we associate clusters of galaxies with dark matter 
haloes we can obtain a(M) from an observational estimation 
of the differential mass function of galaxy clusters, using 
equations @, @ or @, assuming it is well described by 
the PS, SMT or J01 formalisms. In order to obtain a realistic 
a(M) we must ensure that the assumed formalism provides 
an accurate description of the underlying mass function. The 
most reliable results will be obtained for the formalism that 
best describes the mass function, a fact that can be tested 
using numerical simulations. 

3.2 Inverting P(k) from a(M) 
3.2.1 A First try 

As a result of the procedures described in the last section 
we obtain a set of values a(Ri) for i = 1, ...,Nn. For each 
one of them we have an equation that relates this variance 
to the power spectrum A 2 (fc) 



\Ri) = 



A 2 (k)W 2 (kRi)A\n{k) 



(14) 



This is an integral Fredholm equation, with a kernel 
function given by the Fourier transform of the window func- 
tion. With this equation, we are able to recover an estimate 
of A 2 (fc) in a set of values k a a = 1, iVfc. Following Dodel- 
son & Gaztanaga (1999) we will use Greek index for quan- 
tities in fc-space and Latin for those in real-space. 
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Figure 2. Results of the implementation of the A 2 (k e ff) ap- 
proximation, SVD and linear regularization to the functions cr(M) 
obtained from the power spectrum of the tCDM and ACDM mod- 
els 
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-) for the different singular values 



of the matrix W. A great quantity of the ui; are very small and 
are susceptible to roundoff errors. 



A first and rather simple procedure to achieve this is to 
use the fact that equation ( |l4| ) is sharply peaked around a 
characteristic value of k and only a relatively small interval 
of k- values contributes to this integral. Then, a 2 (R) can then 
be understood as an average value of the power spectrum 
in that interval that satisfies the relation ( Peacock, 1991; 
Padilla & Baugh, 2001): 



\R) = A 2 (k eff ) 



(15) 



for some suitable defined value of fc e //. 

In fact, if we assume that the power spectrum is a power 
law A 2 (fc) = Ak n , this relation is exact and fc e // is given by 



ksff = | [9/(n)]i 
where 



(16) 



J(n) 



V n 7 [sm(y) -ycos(y)] 2 Ay. 



It is not necessary to assume that the power spectrum 
behaves like a power law on every scale, as it was mentioned 
earlier. This can be seen by considering that for a given 
value of R only a small range in k is important, and locally 
we can approximate the power spectrum by a power law 
A 2 (fc) oc k n and so equation ( |l6| ) is expected to hold at 
least approximately for this value of n. This is specially true 
in the case of CDM class power spectra, which benefit from 
having slowly varying n(k). 

We can't use this procedure directly to obtain A 2 (fc) be- 
cause equation ( |l6| ) depends on the index n of the unknown 
power spectrum. We can avoid this difficulty by either as- 
suming a fixed value for n on every scale ( e.g. n = — 2 
for scales corresponding to clusters of galaxies) or by ob- 
taining it from a(M). If a 2 (M) behaves like a power law 
a 2 (R) oc R h and 



Figure 2 shows the power spectrum recovered from the 
function a(M) calculated for a rCDM and a ACDM power 
spectrum using this method. The range of masses used in 
this calculation is 10 12 ft _1 M© < M < 10 15 /i _1 M©. It is 
clear from this figure that this procedure gives good results 
in spite of its simplicity. The solutions obtained for A 2 (fc) 
have the correct shape but its amplitude is somewhat un- 
derestimated. 



3.2.2 Singular Value Decomposition 

If we denote by D a and Si the values of A 2 (fc a ) and a 2 (Ri), 
the integral equation ( [Tij ) can be cast as a matrix, yielding 

S =KD (17) 

where K is the N a x Nk kernel matrix and the vectors D 
and S are formed by the values D a and Si respectively. 
Expressed in this way, the first thought about how to obtain 
a solution is by direct inversion of the kernel matrix, that is 



D = K _1 S. 



(18) 



Unfortunately this simple approach does not work. The so- 
lutions obtained inverting K are unstable and often do not 
have a true physical meaning. This due to the fact that we 
try to obtain more information about the power spectrum 
than is contained in a(A'I). This situation causes the matrix 
K to be numerically singular ( although not strictly singu- 
lar). 

In order to obtain the best solution D it is necessary to 
assume that the errors in Si obtained from the mass function 
have a gaussian distribution around their true values Si, 
then the probability to obtain a vector S is proportional to 
exp( 



-^tt) where 



X 



(S-S) T Cg 1 (S-S), 



(19) 
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and Cs is the covariance matrix for S. 

We can relate S to the true value of the power spectrum 
by S =KD. It is easy to see that the solution D that min- 
imizes (|l^) is simply D = K _1 S, and that its covariance 
matrix will be given by Cd = {K j Cg 1 Kj , but, as we 
mentioned before, this solution does not work. To solve this 
problem we use the Singular Value Decomposition technique 
(SVD, see Press et al. 1992). Following Eisenstein & Zaldar- 

riaga (2001), we define K = Cg^K and S =C S ^S, where 

_ I l 
C s 2 is the inverse matrix of Cg constructed by taking the 

square root of the eigenvalues of Cs- We then have 
X 2 = \KD-S\ 2 . 

SVD can be used to find the solution D that minimizes 
this quantity. The singular value decomposition of the K ma- 
trix is given by K — UWV T , where U is a Nr x Nk column 
orthogonal matrix, W is an Nk x Nk diagonal matrix whose 
entrances Wi are called singular values, and V is an Nk x Nk 
orthogonal matrix. With this notation we can express the 
solution (||) as D = VW~ 1 U T S, and the covariance ma- 
trix C p — VW~ 2 V T . If the kernel matrix is singular, then 
some of the elements Wi of the matrix W are zero and the in- 
verse matrix W~ 1 is obviously not defined. But, even when 
none of the singular values are exactly zero, their values can 
be so small as to be dominated by the roundoff errors and 
their inverses w~ 1 will have very large values causing the 
problem to be numerically intractable. Figure 3 shows the 
singular values Wi obtained by performing the SVD of the 
matrix K. It can be seen that there are a large number of 
singular values that are very small causing the problem to 
be singular. 

The solution obtained replacing by zero the quotients 
uij 1 corresponding to the smallest singular values minimizes 
equation ( |is| ) and is therefore the solution D that we are 
looking for ( Press et al., 1992). 

In Figure 2 is also shown the power spectrum recovered 
from a(M) by application of the SVD method for the same 
range of masses than in the u 2 (R) = A 2 (fc e //) approxima- 
tion. It is clear from this figure that SVD gives very good 
results in a wide range of scales, correctly reproducing the 
shape and amplitude of the theoretical power spectrum. 



3.2.3 Linear Regularization 

Another way to obtain a solution for ( |l7[ ) is using a Bayesian 
approach. As we mentioned before the problem we are try- 
ing to solve is singular since we wish to obtain more infor- 
mation about A 2 (fc) from a(M) than that actually avail- 
able. To avoid this lack of information we can supply a 
prior such as the hypothesis of a "smooth" power spec- 
trum. In order to do this we put a second exponential in the 
probability distribution for S, which will be proportional to 
exp \-\ (x 2 + A/3)] , where 



(20) 



Now we have to find the solution D that minimizes the 
weighted sum x 2 +A/3. If the minimization of the functional j3 
alone is non-singular, the general minimization will be so too 
(Press et al., 1992). The first functional ( x 2 ) is a measure 
of the concordance of the solution D with the data S, the 
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Figure 4. Mass function from the GIF tCDM and ACDM sim- 
ulations, determined using FOF halo finder compared with the 
predictions of the PS, SMT and J01 formalisms. 



second term ( 0) can be viewed as a measure of its agreement 
with a prior, implemented by the matrix H, that makes D 
not to vary too much. The A factor in the second term in eq. 
( pp| ) sets the relative weight of the two functionals. The limit 
A — > correspond to the minimization of \ 2 ■ In this work 
we have chosen the matrix H corresponding to a constant 
D (equation 18.5.3 in Press et al., 1992). 

With these definitions the solution for D that minimizes 
X 2 + A/3 can be found by solving the normal equations 

(K T K 



K T S 



XH) D 

and the covariance matrix will be given by 
C D = (^K + XH)' 1 . 

The solutions obtained in this way depend on the value 
of A. This must be determined according to our confidence 
in observational data. A value that gives equal weight to the 
minimization of the functionals x 2 an d /3 is given by 



A = Tr(K T K)/Tr(H). 



(21) 



The results of the application of this method can be 
appreciated in figure 2 where we show the power spectrum 
recovered from the same a(M) used in Kefj and SVD meth- 
ods. The value of A adopted is that of (|2l[). It is clear from 
this figure that Linear Regularization (LR) also provide suit- 
able results for A 2 (fc) in a wide range of scales. 



4 TESTING THE METHOD 
4.1 The GIF simulations 

To test the capability of the method developed in the last 
section to recover the true power spectrum from the differ- 
ential mass function, we applied it to the GIF rCDM and 
ACDM simulations (Kauffman et al., 1999). The values of 
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the cosmological parameters for each simulation are listed 
in Table 1. 

In order to apply our method we first determine the 
mass function of clusters identified in the numerical simu- 
lations. We used a Friends Of Friends (FOF) and a Spher- 
ical Overdensity (SO) halo finders. The mass function de- 
termined in this way is sensitive to the choice of the linking 
length parameter b and density contrast k in the FOF and 
SO algorithm respectively. Both parameters were set in ac- 
cordance with the spherical collapse model as described in 
Eke et al. (1996): b = 0.2 and k = 180 for = 1, and 
6 = 0.164 and « = 324 for Q = 0.3. 

Table 1: values of the cosmological parameters in the GIF 
simulations (Kauffman et al., 1999) 

Model Q.o Oa h <r 8 V Box 

rCDM 1.0 0.0 0.5 0.6 0.21 85ft _1 Mpc 

ACDM 0.3 0.7 0.7 0.9 0.19 140/i _1 Mpc 

The differential mass functions obtained are plotted in 
Figure 4, where we also plot the predictions from the for- 
malisms of PS, SMT and J01 for comparison. As it can be 
seen, the results obtained with the two halo finders are in 
excellent agreement. This allows us to concentrate on only 
one halo finding algorithm which we choose to be the FOF 
algorithm. The later choice is based in the fact that this is 
the algorithm used by J01 in the determination of the val- 
ues given in §2 for the coefficients in the mass function of 
equation (|^) . Another point in support of our choice comes 
from the argument presented by J01 which states that the 
universality of the mass function is more evident using FOF 
than SO. Note that the prediction of the PS formalism fails 
to reproduce the mass functions obtained from the simula- 
tions whereas the SMT and J01 prescriptions give a more 
accurate result. This is no surprise since the GIF simulations 
where used to obtain the parameters in the expressions for 
the mass function in these formalisms. 



4.2 The determination of a(R): 

The next step in the method described in §3 is the deter- 
mination of a(R) from n(M) using equations ( |lo| ) ( p"l| ) and 
(p^[). In doing so, we have only used the mass functions ob- 
tained from the simulations for M > 10 13 h~ 1 MQ, which 
is approximately the range that corresponds to clusters of 
galaxies. Figure 5 shows the results for a(A'I) obtained in 
this way for the rCDM and ACDM simulations, together 
with the a(M) function calculated using (|l4|) and the power 
spectrum present in each model. 

The best solutions for a(M) where obtained by using 
the J01 formalism, indicating that it gives the best descrip- 
tion of the mass function of dark matter halos when these 
are identified using a FOF algorithm. 

The results obtained with the SMT formalism are in 
good agreement with those obtained using the J01 recipe for 
the ACDM simulation, which is not the case with the rCDM 
simulation, showing that this formalism does not describe 
the rCDM mass function at the level of accuracy required 
in our method. The results obtained with the PS formal- 
ism differ notably from those obtained from the theoretical 
model present in the simulations, showing that this repre- 
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Figure 5. Solutions for cr(M) obtained by the application of the 
different formalisms to the mass functions of the GIF simulations. 
For clarity the solutions obtained for the rCDM simulations are 
divided by a factor 1.25. 

sents a poor description of the mass function, a recurrent 
result in the literature on the subject. 

We attempted to determine the covariance matrix Cs 
for a 2 (M) using a bootstrap resampling technique. We de- 
termined n(M) for a large number of samples of groups ran- 
domly selected from the original list from the simulations. 
For each of them we calculate a 2 (M) according to the PS, 
SMT and J01 recipes. Following the notation used in §3.2 
we estimate the covariance matrix by 

(Cfc)y = ((&-&) (Si-Sj)) 

where S represents the theoretical value of a 2 for the model 
considered. 

Figure 6 shows the correlation matrix Corij — 
(C's)ij I \J {Cs)a (Cs)jj obtained with this procedure for the 
J01 solutions determined from the ACDM simulation, the re- 
sults from using the other formalisms are similar. Note that 
all masses smaller (bigger) than Ms are correlated. As the 
mass M only appears in the expressions for a(M) through 
the function G(M) = f™ 8 Mn(M)dM, the value of a(M) 
for masses lesser (bigger) than Ms, have no influence on the 
value of a(M) for masses bigger (lesser) than this limit, so 
the correlations between these two regions should vanish. 
This characteristic is not present in the correlation matrix 
determined by the bootstrap resampling technique, which 
shows spurious correlations that prove the inadequacy of 
this technique for the determination of the covariance ma- 
trix. We have not succeded in determining Cs so we assume 
here a diagonal form. As pointed out by Eisenstein & Zal- 
darriaga (2001) this simplification may cause an underesti- 
mation of the errors in the final power spectrum. For the 
purpose of this work this is a reasonable approximation, but 
it has to be properly taken into account if it is aimed to 
obtain precise constraints upon cosmological models using 
observational data. 
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Figure 6. Correlation matrix of the a 2 (M) function obtained 
from the ACDM mass function using the J01 formalism. This 
matrix was calculated using a bootstrap resampling method. 



4.3 The determination of A 2 (fc) 

Given that in both cases ( rCDM and ACDM) the solution 
ct(M) determined using equation ([L3|) with the general fit 
from J01, is the one that shows a better agreement with 
that of the corresponding cosmological model, we use this 
solution to determine the power spectrum assuming that Cs 
is diagonal. Figure 7 shows the results obtained in this way 
for the rCDM and ACDM simulations using the different 
techniques described in §3.2. For comparison we also plot 
the theoretical power spectrum A 2 (fc) for each model. 

In the implementation of the approximation cr 2 (R) = 
A 2 (k e ff) the value of the index n was not determined from 
a(R) as described in §3.3.1. The solution obtained apply- 
ing this procedure showed an oscillatory behaviour. Instead, 
we fixed n — —2, and on those scales where this approxi- 
mation is valid, the obtained power spectrum shows a good 
agreement with that corresponding to the different models. 
The solution starts to fail where this approximation ceases 
to be valid. The solutions recovered have approximately the 
correct shape, but suffer from the same problem mentioned 
in §3.2.1, underestimating the amplitude of the power spec- 
trum. This analysis shows that this technique could repre- 
sent a powerful tool to determine the power spectrum if the 
index n is determined in a more efficient way. 

The solution obtained applying SVD as described in 
§3.3.2 shows a very good agreement with the correct power 
spectrum over a wider range of scales in the rCDM simula- 
tion but fails to reproduce the correct shape for the ACDM 
case. This behavior was not present in the solution plotted 
in figure 2, where the best results were found using this tech- 
nique. This shows that the solutions obtained using SVD are 
very sensitive to the differences between the used and exact 
a(M), and therefore, to the deviations of n{M) with respect 
to the predictions of the J01 general fit. This indicates that 
the implementation of the SVD method to real data is not 
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Figure 7. Solutions for A 2 (k) obtained with the different tech- 
niques of §3.2, using the solution for <r(M) determined by de 
J01 general fit compared with the correct power spectrum of the 
tCDM and ACDM models. 



very efficient, a situation that may improve with the use of 
the full covariance matrix, instead of a diagonal Cs- 

The best solutions where obtained using the LR tech- 
nique, which reproduces the correct shape and amplitude of 
the power spectrum over a wide range of scales, showing the 
same performance that in §3.2.3. This proves that, at least 
within the hypothesis of a diagonal covariance matrix, this 
is the best technique to recover the power spectrum from 
the mass function of galaxy clusters. 

It should be stressed the fact that the simulations used 
in this test are small in volume. Therefore we have not been 
able to analyze the capacity of our method to deal with the 
presence of rare massive objects. The formalisms of PS, SMT 
and J01 have been extensively tested against the results of 
numerical simulations using volumes much larger than the 
ones used in this work. These analyses have shown that the 
SMT and J01 mass functions account properly for the most 
massive objects. This indicates that the effect of the presence 
of massive clusters in the mass function will be in perfect 
agreement with the theoretical formalisms. Further more, 
the inclusion of these objects in our analysis will improve 
our results and increase the range of scales in which A 2 (fc) is 
recovered, since the numerical problem becomes more easily 
tractable when a wider range of masses is taken into account. 

To test our method in a more realistic situation, we 
also analyzed the effect of introducing errors in the deter- 
minations of individual cluster masses. We constructed 5 
samples of groups from the original one identified in the 
ACDM GIF simulation, by adding gaussian errors to the 
individual masses of amplitudes 10%, 20%, 30%, 40% and 
50%. For each of these new samples we determined the mass 
function. The results obtained show that even an error as 
large as 50% in the individual mass determinations produces 
a small overall effect in the low mass end, causing devia- 
tions of less than 10% in all cases. In the high mass end 
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Figure 8. Solutions for a(M) obtained from the mass function 
of the ACDM GIF simulation by the general fit of J01 fixing the 
values of Q = 0.3, and varying erg in the values 0.8, 0.85, 0.9, 0.95 
and 1.0. Filled circles indicate the solution for erg = 0.9. 
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Figure 9. Solutions for a(M) obtained from the mass function 
of the ACDM GIF simulation using the general fit of J01 fixing 
the values of erg = 0.9, and varying U in the values 0.2, 0.25, 0.3, 
0.35 and 0.4. Filled circles indicate the solution for Q = 0.3. 



( logio(M//i _1 M0) > 14.5) these errors can produce de- 
viations of even 30% due to the small number of objects. 
Nevertheless, their presence does not produce a consider- 
able effect in the final power spectrum, indicating that our 
method is robust even when dealing with large errors in the 
mass estimates. 

4.4 Dependence on and erg 

The method for the determination of the power spectrum de- 
scribed in §3 requires the prior knowledge of the density pa- 
rameter O and the rms mass fluctuations on some reference 
scale, for example R = 8/i _1 Mpc (erg). In this section we an- 
alyze the dependence of our results on these parameters. In 
order to do so we calculate the solution a(M) from the mass 
function of the ACDM simulation (fl = 0.3, erg = 0.9) fixing 
fi = 0.3 and setting erg = 0.8,0.85,0.9,0.95, and 1.0, and 
fixing as = 0.9, and varying ft from 0.2, 0.25, 0.3 , 0.35 and 
0.4. The results obtained in this way are plotted in Figures 
8 and 9. 

Figure 8 shows that except for the high mass end, the 
effect of changing erg in the solutions of a(M) is simply a 
change in the overall normalization. On the other hand, the 
shape of the solutions is affected in the high mass limit. As is 
shown below, this effect propagates into the power spectrum 
causing differences in the low k end. 

Figure 9 shows that the effect of changing the value 
of Q can be more serious, as the solutions obtained have 
different behaviors than that observed in the £1 = 0.3 case. 
This will cause differences in the obtained power spectrum 
in all scales. 

In order to analyze the effect of these variations on the 
power spectrum we use the LR technique to obtain A 2 (fc) 
from the cr(M) functions plotted in Figures 8 and 9. In doing 
so we used the same value of A used in §4.3 when obtain- 
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Figure 10. Solutions for A 2 (fc) obtained aplying the LR tech- 
nique to the solutions for cr(R) from the A CDM simulation by 
fixing Q and varying erg (see Figure 8). 



ing the power spectrum from this simulation. The results 
obtained this way are plotted in Figures 10 and 11. 

Figure 10 shows that, as we pointed out previously, the 
effect arising from changing the value of erg is important 
in the lower values of k, where it can produce considerable 
changes in the solutions. For the higher values of k, the 
solutions obtained are not very different for the one obtained 
for the correct value of erg. This is the same effect observed in 
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Figure 11. Solutions for A (fc) obtained aplying the LR tech- 
nique to the solutions for a(R) from the A CDM simulation by 
fixing erg and varying S7 (see Figure 9). 
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Figure 12. Results obtained for the power spectrum P(k) of 
the A CDM GIF simulation ussing the LR technique and the 
FKP method to the same group sample, with logarithmic (upper 
panel), and linear (lower panel) scales. 



the analysis of a(M) where we only found a renormalization 
and not a change in shape. 

The most important effect on the power spectrum, is 
that caused by the use of a wrong value of fi. As it can 
be seen in Figure 11, the solutions obtained differ notably 
between them, and show different behavior in every scale. 
The reason for this is that the effect of changing Q in a(M) 
is not a simple renormalization, but a change in the slope of 
the solutions. This makes f2 the most important parameter 
in the model. 
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5 COMPARISON TO THE FKP METHOD 

As a final test of our method we compare the outcomes of 
the method presented in this work, and that of FKP, as 
described in Zandivarez, Abadi & Lambas (2001). 

We applied both methods to the GIF simulations and 
the results are plotted in Figures 12 and 13 corresponding to 
the ACDM and rCDM models respectively, where we plot 
the power spectrum P(k), output of the FKP technique, in 
logarithmic (upper panel) and linear scales (lower panel). In 
order to compare these results we have arbitrarily normal- 
ized the power spectrum derived from the FKP technique 
to that of the linear mass power spectrum, which is the out- 
come of our method. 

By inspection to these figures, it is clear that FKP and 
our method yield similar results, showing the power and per- 
formance of our approach. We notice that the solutions for 
P(k) obtained from the mass function show a much better 
agreement with the power spectrum of the corresponding 
model than those obtained using the FKP technique. The 
range of validity of our solutions must be kept in mind, since 
the apparent good description of P(K) in the higher values 
of k correspond to the region where LR gives a constant 
solution for A 2 (fc), which in turn gives P(k)ak~ 3 , the ap- 
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Figure 13. The same as Figure 12, but for the r CDM GIF 
simulation. 

proximate behavior of the power spectrum in these scales. 
Even though LR gives a more accurate description of P(k) 
in the range where it has been proven to be applicable, the 
estimated power spectrum using the FKP technique works 
for higher values of k. As the range of validity of our so- 
lution depends on the mass range of the differential mass 
function, the inclusion of low clusters masses will improve 
this situation. 

As it can be seen in the lower panels of Figures 12 and 
13, the solution obtained by the use of LR gives a very accu- 
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rate description of the power spectrum for the lower values 
of k, where the FKP method fails to give a satisfactory an- 
swer. Besides, the logarithmic bins in LR solution's sample 
these values of k in a better way than the linear bins of 
the FKP technique, showing that our method is better in 
studying the power spectrum at these scales. 

It is important to stress the fact that our method pro- 
vides the linear power spectrum, while the solution of the 
FKP technique is affected by nonlinear effects. These are not 
important in the range of scales analyzed here, but this may 
be significant when objects of lower masses are included in 
the analysis. 

A significant and important difference in our technique 
is that it naturally allows the use of information from differ- 
ent scales and objects, ranging from poor to rich clusters of 
galaxies in an homogeneous way. On the other hand the FKP 
technique will give different results for subsets of systems of 
galaxies with diferent masses, reflecting their relative bias 
with respect to the underlying mass distribution. Further- 
more, in our method it is desirable the inclusion of objects 
with a very broad range of masses, since this allows for a 
better determination of the power spectrum over a wider 
range of scales. 



6 CONCLUSIONS 

In this work we have developed a new method to determine 
the linear mass power spectrum from the mass function of 
galaxy clusters based on the PS, SMT and J01 formalisms. 
The rms fluctuation of mass a(M) is determined using the 
expression for the mass function in the form of a differen- 
tial equation. We then use different techniques from integral 
equations theory to obtain the dimensionless power spec- 
trum A 2 (fc) from o(M). 

The implementation of this method to the GIF rCDM 
and ACDM simulations shows a high efficiency in recovering 
the correct shape and amplitude of the power spectrum in 
a wide range of scales, showing a high performance and real 
practical utility. Although the different solutions for a(M), 
obtained from the simulations, show an overall agreement, 
it is clear that the J01 prescription gives a much better re- 
sult than the PS or SMT formalism, showing that the J01 
general fit gives an excellent description of the differential 
mass function of clusters identified in the simulations. 

In the theoretical case analyzed in §3.2 the best solu- 
tion was that obtained from the implementation of SVD, 
which recovers the power spectrum in a wide range of scales. 
However, from the analysis of the results obtained from the 
simulations we conclude that the best technique to invert 
the integral equation in our problem is linear regularization. 
This situation can change with the use of the full covariance 
matrix and so, it must be analyzed without the assumption 
of a diagonal Cs- 

The comparison of the results obtained with our method 
with those obtained by the standard FKP technique, shows 
a very good general agreement. Nevertheless, in the scales 
of validity of our solutions, our answer gives a much more 
accurate description of the power spectrum than the ones 
obtained by the use of the FKP method. This is specially 
valid in the lower values of k where the FKP solution fails 
to give the correct answer. Further more, the logarithmic 



bins of our method sample these scales in a more efficient 
way than linear bins (as is the case in the FKP output). This 
produces important effects in the analysis of the power spec- 
trum in large scales. The FKP technique works for higher 
values of k than the LR solution. However, this can improve 
significantly if information of n(M) for lower masses is in- 
cluded in the analysis. 

In order to apply our procedure to real data we need 
an estimate of the mass function of galaxy clusters. Un- 
fortunately, present estimates of this important statistical 
measure are not completely satisfactory since quoted values 
are significantly different and moreover they are referred to 
different cluster mass definitions (see for instance Bahcall & 
Cen, 1993; Biviano et al., 1993; Girardi et al., 1998; Girardi 
& Giuricin, 2000; Reiprich & Bohringer, 2001; Bahcall et al., 
2002). 

As we have already pointed out, even in the case where 
errors in the individual masses are large, the overall effect in 
the mass function is not important when compared to other 
sources of uncertainties such as the use of different mass def- 
initions. Different methods for cluster mass estimates ( i.e. 
optical, x-ray, gravitational lensing) are clearly correlated, 
but provide different values for the mass of the same clus- 
ter. Furthermore, it is not clear how these mass definitions 
are related to the ones used in numerical simulations ( i.e. 
FOF b = 0.2, FOF b = 0.17, SO n = 200, etc.) which also 
show differences among them ( see White, 2001). This is of 
extreme importance and is the principal source of uncertain- 
ties in our method. Given that observational cluster masses 
are estimated using one of the methods mentioned above, it 
may be incorrect to assume a fit to the mass function whose 
coefficients have been determined using FOF masses, such 
as the J01 formalism. 

This indicates that the use of expressions derived from 
theoretical formalisms in the analysis of observational data 
( a common practice in the literature) must not be made 
without a careful analysis of the relationship between the 
different mass definitions involved. 

There are other statistics available in the literature, 
such as the temperature function (Henry & Arnaud, 1991, 
Markevitch, 1998 , Blanchard et al., 2000; Ikebe et al, 2001) 
of X-ray samples. The high correlation between tempera- 
ture and mass suggested by numerical simulations (Evrard, 
1990; Evrard, Metzler & Navarro, 1996; Yoshikawa & Suto, 
2000, Mathiesen & Evrard, 2001), has lead many authors to 
use temperature functions in confronting observations with 
theoretical predictions, instead of mass functions. 

The method described here can equally be applied to 
the temperature function of galaxy clusters, but, in doing 
so, we need a reliable mass-temperature relation. Unfortu- 
nately this transformation has not been yet determined in 
an accurate way, and different T — M relations give equally 
different results between them. 

The observational determination of the mass function of 
galaxy clusters can improve in the next few years with the 
outcome of new mass estimates from gravitational lenses, or 
even from virial mass estimates obtained from the new large 
catalogs that will be soon available such as 2dF and SDSS. 

If the above discussion is properly addressed the meth- 
ods described in this work can be very useful tools to provide 
new and reliable determinations of the linear mass power 
spectrum. 
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APPENDIX A: OBTAINING a{M) FROM THE 
PS, SMT AND J01 FORMALISMS 

The expression for the mass function in the PS, SMT and 
J01 formalisms can be re-interpreted as differential equa- 
tions to obtain o(M ) from n(M). If we use the PS formalism 
then the mass function is given by: 



n(M, z)dM : 



7T M 2 <j(M,z) 
that can be written as 



dln((j) 



dln(M) 



2cr 2 (M,z) 



dM 



Mn(M)dM 



[l_t e -{St)da dM 
V it a 2 dM 



where all dependence on a is in the right hand side. Integrat- 
ing in both sides with respect to M between a mass M and 
Ms — ^p(8h~ Mpc) corresponding to a sphere of uniform 
density and radius R — 8/i _1 Mpc, we get 



G(M) 



where p 



M 8 



Mn(M)dM : 



and (i$ 



dp 



(Al) 



Using the definition of the error function 
erf(a;) = — = / dp 



and using the value of as as a constant of integration, (Al) 
can be expressed as 



G(M) = p 



erf(- 



4 



erf(- 



"^j2o % V2a(M)' 

and then we obtain the expression for <r(M) given in equa- 
tion @ 

a(M) = Sc 



erf( 



V2a 



G(M) 



If we assume the validity of the SMT formalism then 
the differential equation for a(M) is 



Mn(M)dM = 



ty a 1 



1 + 



V \fa5 c 



2p 



"(f^) da_ 



dM 



dM 



This is a little more complicated than the corresponding 
to the PS formalism but can be solved in an analogous way. 
Integrating between M and Mg we get 



G(M) = Ap— 



e-» dp + 2~ v 



-2p - 

p e 



dp 



(A2) 



where G(M) is given by (Al) and now p = \/~% fjjry 

The first term of the right hand side in (A2) is the 



same than in the PS formalism and is equal to erf(zj) , 
the difference in this case is given by the second term. If we 
use the definition of the incomplete gamma function 



P(a,x) = 



r( ) 



L dt 



a > 
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we can see easily that 



-2p - 

p *e 



dp 



= --p p --p.mm 



e 



1 



Then, if we define the function 

*(«) = erf (x) + Q - p) P (i - p, : 

we can write ( |A2| ) as 

5 C N , / 5 C 



G(M) = 



%/2. 



<T8 



")-*( 



from where we obtain the expression for cr(M) in the SMT 
formalism given in equation ( |ll| ) 

a(M) = 



V2cr 8 ' 



pA 



It is possible to consider that the J01 mass function 
gives a more accurate description of the real behavior of the 
mass function. If this is the case, then we use the expression 

which must be re-interpreted as a differential equation for 
a(M) as a function of n(M, z). This expression can be writ- 
ten as 

Mn(M)dM = Ape- i[s+Bn —dM 
v ' ' dM 

where S = ln(<r -1 ) = — ln(cr). Integrating in both sides with 
respect to M between M and M% and keeping the definition 
of G(M), we get 

G{M)=Ap(p{Ss)-*(S)) 
where 



*(S) = 



\x+B\' 



dx 



Then we can write 



VP" 1 



*(5 , 8 



G(M) 
Ap 



from where we get the following expression for a(M) 

G(M)~ 



a(M) = exp j-* 1 
which is given in equation (|l3|). 



Ap 
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